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Abstract 

We formalize the notion of sampling a function using 
k-d darts. A k-d dart is a set of independent, mutually 
orthogonal, /c-dimensional subspaces called k-d fiats. 
Each dart has d choose k fiats, aligned with the co- 
ordinate axes for efficiency. We show that k-d darts 
are useful for exploring a function's properties, such as 
estimating its integral, or finding an exemplar above a 
threshold. We describe a recipe for converting an al- 
gorithm from point sampling to k-d dart sampling, as- 
suming the function can be evaluated along a k-d fiat. 

We demonstrate that k-d darts are more efficient 
than point- wise samples in high dimensions, depend- 
ing on the characteristics of the sampling domain: e.g. 
the subregion of interest has small volume and eval- 
uating the function along a flat is not too expensive. 
We present three concrete applications using line darts 
(1-d darts): relaxed maximal Poisson-disk sampling, 
high-quality rasterization of depth-of- field blur, and es- 
timation of the probability of failure from a response 
surface for uncertainty quantification. In these appli- 
cations, line darts achieve the same fidelity output as 
point darts in less time. We also demonstrate the ac- 
curacy of higher dimensional darts for a volume esti- 
mation problem. For Poisson-disk sampling, we use 
significantly less memory, enabling the generation of 
larger point clouds in higher dimensions. 

1 Introduction 

In many applications we are interested in estimating 
some global property of a function because it is diffi- 
cult to calculate that property exactly. Sampling is the 
process of randomly selecting samples^ subsets of a do- 
main. The function is evaluated at these subsets, and 
the global property is estimated based on those values. 

In typical sampling processes, the samples are points. 
However, a recurring challenge is to deal efficiently with 
the case that the interesting part of the domain is very 
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Figure 1: Sampling long and thin subregions (gray) 
using points (left), lines (center), and planes (right). 
Point samples may be cheap to generate and evaluate, 
but they contribute nothing to the final result if they 
miss the region of interest. Misses (blue) are frequent 
for regions with a small volume. Samples of higher di- 
mensions, or k-d darts, often intersect (red) the region 
of interest, especially if the region is long and thin. A 
k-d dart's greater expense is offset by it providing more 
information. 

small compared to the entire domain. For example, 
suppose we have a function over a domain, and we are 
interested in estimating the volume of the subdomain 
where the function is negative. If this subdomain has 
a very small volume, only a correspondingly very small 
fraction of uniform sample points will land in it; see 
Figure 1. Consequently, point sampling will require a 
large number of samples to get any estimate, and will 
be inefficient since most samples will not contribute. 

We propose the k-d dart to address this problem. 
One key idea is that rather than evaluate the function 
at a single point, we evaluate it in a higher-dimensional 
region. For each sample, we evaluate the function along 
a set of higher-dimensional flats (i.e. lines, planes . . . 
hyperplanes). The second key idea is to use a set of 
mutually orthogonal flats, aligned with the coordinate 
axes; a rC- d dart denotes this set of flats. Randomly 
oriented flats have been considered before, but orthog- 
onal flats are more efficient and have better worst-case 
performance when probing high aspect-ratio settings. 
To ensure that the expected mean of the function esti- 
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mates is correct, each of these fiats is chosen indepen- 
dently. An important case of flats are one-dimensional 
lines. Using our previous example, we may find the 
points along the line where the function value is zero, 
then partition the line into segments where the func- 
tion value / is strictly positive or negative, and finally 
estimate the volume where / < from the negative- 
interval lengths. While these samples are more expen- 
sive to compute, they are more powerful; depending on 
the function they can generate better results for the 
same amount of effort. 
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(a) Monte Carlo Sampling (MC) 
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(b) Latin Hypercube Sampling (LHS) 

Figure 2: Estimating the volume of a sphere using ran- 
dom sampling via k-d fiats, k = 0, 1, 2. For each sam- 
ple size, we performed 100 experiments and calculated 
the RMS error. The reported CPU time is the total 
time consumed by these experiments. For MC sampling 
(a) plane samples consumed an order of magnitude less 
time to achieve the same error as point samples. The 
savings were even more for LHS (b). 



A simple example that demonstrates this concept is 
estimating the volume of a unit sphere by sampling 
from its bounding box: / = — 1 inside the sphere and 



outside it, and we seek an estimate of J^^q 1. Figure 2 
shows the relation between error and time as the size 
of the sample increases. We sampled using k-d fiats of 
dimension /c = 0, 1, 2. For a point sample, we checked if 
the point was inside the sphere. For higher dimensions, 
we calculated the fraction of the fiat inside the sphere; 
see Figure 3. We performed both Monte Carlo (MC) 
and Latin Hypercube Sampling (LHS). For each sample 
size we ran 100 experiments and calculated the error 
in the volume estimate. Plane samples consumed less 
CPU time than point samples for the same RMS error. 
For MC sampling the payoff was about a factor of 5, 
and for LHS sampling the payoff was 3 to 8 orders of 
magnitude! The reasons behind these gains are the 
following: 

• Evaluating / along k-d fiats is cheap; in this case 
we exploited the analytic function of the sphere. 

• A k-d fiat gives more information as k increases. 

• A fiat is cheap to generate. Each k-d fiat requires 
d — k random numbers; here d = 3. 



(a) 0-d flats 



(b) 1-d flat 



(c) 2-d flat 



Figure 3: k-d fiats used to estimate the volume of a 
sphere. The fraction of the sub-fiats inside the sphere 
estimates the function average. 



In general, evaluating the integration function along 
a k-d fiat costs more than at a single point. However, 
for many problems, this extra cost is offset by the supe- 
rior capability of a k-d fiat to capture narrow regions. 
For instance, consider Figure 4(a), where a line fiat per- 
pendicular to that narrow region of interest will cap- 
ture it regardless of its thickness. On the other hand, 
the probability of a point sample landing in the region 
approaches zero as the thickness decreases. 

The purpose of this paper is to formalize and demon- 
strate the k-d dart approach. In Section 4, we show 
three practical applications: completing a relaxed max- 
imal Poisson-disk sampling, in dimensions 4-30; ren- 
dering depth-of-field blur in four dimensions; and esti- 
mating the probability of failure from a response sur- 
face for uncertainty quantification, where the proba- 
bility is small, e.g. le-5, and the dimension is large. 
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e.g. 15. In each of these three apphcations the space 
is of moderate dimension, and hne darts are particu- 
larly effective. The experiments in Section 5 verify the 
accuracy of higher-dimensional darts, using a volume 
estimation application. 

2 Previous Work 

Our work generalizes sampling. There are many pat- 
terns for generating samples. For graphics, maximal 
Poisson-disk sampling (Section 2.1) is common. Much 
prior work focuses on dimensions 2 and 3, but our inter- 
ests extend beyond. Rendering applications use sam- 
pling in many forms (Section 2.3), often with high di- 
mension. The field of uncertainty quantification (non- 
graphics. Section 2.4) attempts to quantify the range 
of a computation, typically by performing the compu- 
tation many times with a well- distributed selection of 
input parameters. This is important in computational 
science for predictions and reliability. 

Our k-d dart is a particular conception of high di- 
mensional sampling. Line sampling has been explored 
in varied contexts already, within graphics and other 
domains. Unfortunately much of this work is isolated 
and does not consider dimension as a free parameter; 
we hope to help provide a unified view of these ap- 
proaches. 

Lines were used early in the history of Monte Carlo 
sampling. For instance, "Buffon's needle problem" was 
published in 1777. It considers the number of inter- 
sections of randomly-oriented short line segments with 
axis- aligned regularly spaced infinite lines. Its solution 
involves geometric probability, either through integral 
geometry [25] or MC experiments [12], and gives an es- 
timate of the value of tt. Both the needle problem and 
volume estimation by line darts use finite objects, vary- 
ing orientations, and infinite probes; but which things 
are known, measured, and estimated are different, as 
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Figure 4: Extreme subregion shapes. In general, a 
k-d fiat has a better chance to intersect a region of in- 
terest as k increases. For a given region volume, the ad- 
vantage is higher for stretched regions than for square 
ones. 



are which quantities are uniform-random and which 
are uniform-deterministic. For example, using sets of 
orthogonal fiats means the indices of the fixed dimen- 
sions of our fiats are evenly spaced deterministically, 
whereas in Buffon's problem the geometry of the rule 
lines are evenly spaced. 

In neutron transport physics simulations, a class of 
Monte Carlo algorithms known as "track length esti- 
mators" essentially performs Monte Carlo estimation 
using line segments [28] . This is in contrast to the "col- 
lision estimators" class that estimates using point sam- 
ples. In graphics, these collision estimators correspond 
to standard volumetric photon mapping, estimation 
using photon scattering locations. Also track-length 
estimators (using line samples) correspond to photon 
beams; estimation uses entire random- walk path seg- 
ments [14]. In surface reconstruction and CAD mod- 
eling, we can count the intersections of unoriented line 
samples with the surface, then make use of the integral 
geometry Cauchy-Crofton formula to estimate integra- 
tion quantities such as surface area or enclosed vol- 
ume [19]. To get the right estimate with low variance, 
it is crucial to select the sample lines using the right 
probability model. For example, models for bundles 
of uniformly-spaced parallel lines, reminiscent of both 
k-d darts and the regularly spaced lines in Buffon's 
needle problem; models for chordal lines; and pseudo- 
random and other numerical sequences have all been 
studied in the context of sampling surfaces [26]. 



2.1 Relaxed Maximal 
Sampling 



Poisson-Disk 



Maximal Poisson-disk Sampling (MPS) is a popular 
graphics technique to distribute a set of points in a 
domain. The points are random and have a blue-noise 
spectrum, which is well-suited to the human visual sys- 
tem and helps avoid visual artifacts. The points have 
a minimum distance between them, r/, which helps to 
use the point budget efficiently. We denote the maxi- 
mum distance from a domain point to its nearest sam- 
ple by Tc. In a maximal sample, the Vf disks around 
the points overlap to cover the whole domain, leaving 
no room to add another point, and Tc < r/. Otherwise, 
maximality is relaxed, and we measure how far the ge- 
ometry is from maximality by the distribution aspect 
ratio €r = Tc/vf > 1. Point clouds with a meaningful 
upper bound on are separated-yet-dense, also known 
as "well-spaced." 

MPS algorithms abound, and often achieving maxi- 
mality is the most challenging part. In some applica- 
tions maximality is not required [17]. The acceptable 
relaxation of maximality depends on the application. 
For example, in Voronoi mesh generation [4], the cells 
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have an aspect ratio bound that varies smoothly with 
the relaxation, 2e^. Some methods sacrifice maximal- 
ity to terminate more quickly, but explicit statements 
about the achieved are rare. 

Many methods use some form of a background grid 
for point location and proximity queries [6, 15, 32]. The 
background grid may be refined, as in a quadtree, to 
track the remaining voids (uncovered regions) [7, 33]. 
This can be made efficient in dimensions up to about 
5 [5]. However, even in dimensions below 5, refinement 
methods can run out of memory as the sample size in- 
creases. Memory problems are exacerbated on a GPU. 

Uncertainty quantification motivates MPS sampling 
in higher dimensions, e.g. 10-30. No MPS methods in 
the literature scale to these dimensions due to the so- 
called "curse of dimensionality." Classical dart throw- 
ing [2, 3] is not strongly dependent on dimension, but 
as the number of accepted samples increases, the run- 
time for the next sample becomes prohibitive and the 
algorithm must terminate well before maximality. 

Consider sampling a unit box with disk radius r in 
dimension d. Some issues are fundamental to the prob- 
lem, independent of any specific algorithm or applica- 
tion. The size of a maximal sampling n is indetermi- 
nate, but its lower and upper bounds grow exponen- 
tially in 1/r and d. The kissing number, the number of 
disks that can touch another disk, grows exponentially 
in d. (The geometry literature discusses these issues 
extensively. The densest packings and largest kissing 
numbers by dimension are summarized in "A Catalogue 
of Lattices" [23].) 

The curses of dimensionality for grid-based methods 
go beyond those unavoidable issues: 

1. The base grid size grows exponentially and faster 
than the output point set. 

2. A grid cell is refined into 2^ subcells. 

3. The number of nearby cells that might contain a 
confiicting sample grows faster than the kissing 
number. 

4. The ratio of misses/hits grows exponentially with 
the dimension [5]. 

Item 1 arises because the size of the base grid is usu- 
ally chosen such that each cell can accommodate at 
most one point: base squares have diagonal length r 
and edge length r/y/d. Worse, at maximality, the num- 
ber of empty base cells grows exponentially with the 
dimension. These empty cells increase the time, and 
especially the memory requirements, to prohibitive lev- 
els. One possible solution is to choose a base grid level 
with edge length 2d^ so that every square must contain 
one point. However, due to Item 2, this approach just 



defers the problem until cells are required to be refined 
a couple of times to represent voids. Because of the 
kissing number, representing voids through geometric 
constructions [6] or explicit arrival times [15] do not 
appear to be viable solutions at high dimension. 

Some approximate methods [32] put an upper bound 
on the number of misses per cell, which partly addresses 
Item 4. The drawback is that the sampling is not max- 
imal. How far it is from maximality has not been an- 
alyzed, but volume arguments suggest that for a fixed 
box size, the number of allowed misses must grow ex- 
ponentially in d to bound the linear distance between 
an uncovered point and a sample's disk. 

While some of these issues affect runtime, the real 
curse is the memory requirements for quadtrees. Sim- 
ple MPS [5] is the quadtree method with the best mem- 
ory scaling by dimension. It seems unlikely to extend 
to even = 10 in the near future. Simple MPS uses a 
fiat quadtree of same-sized squares, periodically refin- 
ing all remaining squares uniformly. We compare our 
results to Simple MPS in Section 4.1.4. 

2.2 Motivation for High-Dimensional 
Point Clouds 

In the design of computer experiments, we generate 
points in parameter space, then evaluate a function at 
the sample points. We often build a surrogate model 
based on those points' values, e.g. Kriging models. The 
dimension is equal to the number of parameters, so 
can be very high. The time it takes to generate the 
points is often very small compared to evaluating the 
function, so it is worthwhile to spend the time to find 
a set of points that span the space efficiently. Well- 
spaced points use the point budget efliciently, and pro- 
vide bounds on the condition number and interpolation 
error. 

For some applications, if the function is not too ex- 
pensive, it may make sense to sample the function di- 
rectly using k-d darts. For example, if the integral 
of a function over a ffat is available analytically, then 
sampling it directly using k-d darts would be very ef- 
ficient. Otherwise, k-d darts can provide well-spaced 
points (see Section 4.1); we can build a surrogate model 
over those points; then k-d darts can integrate the sur- 
rogate model analytically. 

Another application where sample points are re- 
quired is finite element simulation, where we need a 
computational mesh of the points. 

2.3 Rendering 

High-quality rendering is an important application of 
multi-dimensional sampling. Photorealistic effects like 
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motion blur, depth-of-field (defocus) and soft shadows 
can be expressed as integrals over multiple dimensions. 
Classical techniques often employ stochastic point sam- 
pling to estimate these integrals. Noise-free rendering 
using point sampling can require a large number of sam- 
ples, which can be extremely expensive for complicated 
scenes. 

Thus a long history of research has targeted choos- 
ing samples wisely. Mitchell's classic antialiasing pa- 
per [22], for instance, "focus [es] on reducing sampling 
density while still producing an image of high quality," 
primarily in the context of ray tracing. Metropolis light 
transport [31] "performs especially well on problems 
that are usually considered difficult" by using muta- 
tions to preferentially (but in an unbiased way) sample 
light paths in interesting regions of the path space. 

Besides choosing samples carefully, another approach 
toward the same goal is to reduce the number of re- 
quired samples through techniques such as sample reuse 
and/ or caching. For example, a notable recent advance 
in this area, by Lehtinen et al. [18], specifically notes "a 
clear need for methods that maximize the image qual- 
ity obtainable from a given set of samples" and exploits 
anisotropy in the temporal light field to reuse samples 
between pixels. Our work has a similar goal — making 
the best use of a limited number of samples — but in- 
stead reduces the number of necessary samples by in- 
creasing their dimensionality. 

k-d darts formalize a general way of multi- 
dimensional sampling. Several early graphics applica- 
tions use multi-dimensional samples in limited and spe- 
cific scenarios. The OpenGL accumulation buffer [11] 
uses a form of k-d darts for motion blur: each out- 
put pixel is an aggregate of several input pixels, each 
of which may span a 2-d region (pixel area) for a con- 
stant shutter time. Essentially, the accumulation buffer 
samples 2-d x-y images in a 3-d x-y-time space. Nelson 
Max [20] uses a scan-line visible surface algorithm that 
generates line samples, describing how to use the infor- 
mation in these samples to create antialiased images. 

Recent research has also shown promise in rendering 
high-quality motion blur using multi-dimensional sam- 
ples. Gribel et al. [10, 9] present the use of hne samples 
(our "1-d fiats"). In their 3-d domain (x,7/,t) they fix 
X, y and perform a 1-d fiat in the t domain^. They also 
extend their implementation to render motion-correct 
ambient occlusion. 

More recently, line samples have proven useful in the 
representation of light. Sun et al. [29] represent light- 
ing and viewing rays directly in a 6-d Pliicker space, 
which allows an efficient formulation of finding nearby 
lighting rays. This, in turn, allows accurate, fast ren- 

-"^In their 2011 paper, Gribel et al. use 2-d darts but call them 
line samples because they are lines in x-y space. 



dering of large scenes with single scattering even in the 
presence of occlusions and specular bounces. The pre- 
viously mentioned work of Jarosz et al. [14] concen- 
trates on better representations for light paths in pho- 
ton tracing for the purposes of rendering participating 
media in light interaction; previous methods had used 
photon particles. Instead, they represent and store full 
light paths (samples) , resulting in a more compact and 
expressive lighting representation with corresponding 
performance benefits. 

Jones and Perry [16] experimented with using ana- 
lytical line sampling for anti-aliased polygon rendering. 
They shoot single-dimensional darts across a pixel's 
surface, analytically compute triangle coverage for each 
of them, and then average them to obtain pixel colors. 

Our paper generalizes the idea of multi-dimensional 
darts for sampling, and it is this generalization that 
helps us design renderers that use A:-d darts in differ- 
ent configurations. For rendering depth-of-field effects 
with added anti-aliasing, we present a generalized con- 
figuration using a full 1-d dart in Section 4.2. That 
is, we use a set of orthogonal fiats rather than just a 
single flat as in prior work. We compute high-quality 
depth-of-field images efficiently. 

While we demonstrate just one configuration of 
k-d darts, several different strategies are possible for 
depth-of-field and other effects. In general, k-d darts 
offer a sampling process that converges faster and has 
lower noise than point sampling; the use of k-d darts of- 
fers the potential benefit of faster convergence but must 
be weighed against the higher complexity per k-d dart. 

2.4 Monte Carlo Sampling 

Random sampling is one of the oldest [12] and most ro- 
bust methods for uncertainty quantification. The prin- 
cipal use of Monte Carlo (MC) sampling is to approxi- 
mate a high- dimensional integral with a sample mean. 
The primary drawback of MC is the slow rate at which 
the sample mean converges to its true value. The stan- 
dard error in the computed mean for n samples is 

C^err = —/= • (1) 

Although this rate of convergence in n is very slow, 
the number of dimensions, does not appear in Equa- 
tion (1). MC's primary advantage is that it is not sub- 
ject to the curse of dimensionality. Significant effort 
has been invested in developing variants of MC with 
faster rates of convergence; a full review is out of scope. 
Latin Hypercube Sampling (LHS) [21, 24] is known as 
"N-rooks sampling" in the graphics community. Im- 
portance sampling [8] involves preferentially sampling 
rare cases and compensating by reducing their weight. 
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3 Dart Framework 

Given the ideas of the method (Section 1) and where it 
might be useful (Section 2), we now define it formahy 
and give a general recipe for converting a point sam- 
pling algorithm to a k-d dart sampling algorithm. We 
consider two application scenarios. The first is sam- 
pling to evaluate a function, as in numerical integra- 
tion; see Section 3.1. The second is sampling to find 
locations in the domain where the function has a par- 
ticular value; see Section 3.2. 

We begin with terminology. Function / is defined 
over domain X C M^. Often X = [0, 1]^. A flat F is a 
A:- dimensional subspace of the domain, the vector space 
defined by fixing d — k coordinates. A k-d dart is a 
sample defined by a set of / = (^) /c-dimensional fiats, 
s^, one for each combination of fixed coordinates. So 
in M^, a 1-d dart might consist of the two 1-d fiats at 
X = 0.5 and y = 0.5. 

3.1 Evaluating a Function over a Do- 
main 

Consider an algorithm that receives a single value, 
y = /(x), from evaluating a function / at a uniform 
random point x. To get an analogous single value y^ 
from a k-d dart, we calculate the average of / over its 
fiats, weighted by the relative measure of the fiats. We 
generate a k-d dart's fiats one at a time. For each fiat 
F we choose its fixed coordinates uniformly at random, 
then: 

1 . Clip F at the boundary of the domain and estimate 
its relative volume |i^| = /pi- For the unit box 
this is trivial. 

2. Integrate the function along the clipped fiat: G = 

3. Now the weighted average of / over the fiat is 
simply the ratio of these two computations: H = 
G/\F\. 

Then y^ is simply the average of H over the flats: 
y' = Y^l^i Hi/ 1. We can combine multiple k-d dart y' 
measurements in the same way that we might combine 
measurements made at multiple point samples. Thus 
we can convert point-sampling function evaluation to 
k-d dart function evaluation. 

A generic approach (for any /c-dimensional flat) is to 
use numerical integration over a discretization of the 
flat. Represent a flat with a uniform grid (mesh). (We 
assume it is too expensive to mesh the entire domain.) 
In (1), clip the grid elements at the domain boundary, 
generating simplices. In (2), evaluate / at the grid 



points and perform standard numerical integration over 
the grid simplices. 

Our depth of field, probability of failure, and volume 
estimation applications follow this recipe, because they 
use Monte Carlo integration, but the integration along 
the line fiats can be made faster and more accurate 
given application-specific knowledge about the func- 
tion. For probability of failure, the function we inte- 
grate is the 0-1 indicator function of failure, rather than 
the continuous value of the response function. We im- 
prove accuracy by finding the roots of a single- variable 
equation to find the boundary points where the indi- 
cator function switches its value. For depth of field, 
/ is the contribution of one ray (photon) to a pixel, 
and we seek to estimate the contribution of all photons 
over all focal depths for each pixel. For efficiency we 
use discrete algorithms to find the occlusion boundaries 
of triangles, then integrate a continuous function over 
each non-occluded segment of a triangle. 

3.2 Identifying Points in the Domain 
with a Particular Function Value 

MPS represents a diff'erent category of application. In 
this application, instead of integrating /(x), we are 
looking for a random x that satisfies f{x) > 0, i.e. 
finding a point outside all prior disks. Point-sampling 
techniques would sample the domain until such a point 
was found, which is expensive if the region of interest 
is small compared to the size of the domain. Instead, 
we can use k-d darts for sampling using the following 
generic recipe: 

1. Clip F at the boundary of the domain. Retain 

a representation of the clipped fiat. (For MPS, we 
clip a line by the unit cube.) 

2. Retain the portions of g where the function has 
the particular value of interest, e.g. / > 0. (For 
MPS, these are the subsegments outside all prior 
disks.) 

3. Return a random point from g. 

4 Applying the Framework 

We now describe how to apply the k-d dart framework 
to our representative applications in more detail. 

4.1 Relaxed MPS 

Our relaxed maximal Poisson-disk sampling algo- 
rithm is a variant of traditional dart-throwing, throw- 
ing point darts and keeping those that hit an uncov- 
ered region (void). Algorithm 1 specifies our imple- 
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Algorithm 1 A classical dart throwing algorithm us- 
ing line darts. 

while maximality estimates are inadequate do 
generate a line dart 
for all i =permutation(l..(i) do 

generate line segments g = sjo domain 
for all samples p do 

subdivide g = g\D{p) 
end for 
if ^ 7^ then 
count s-"^ as a hit 

select a sample point uniformly from g 
skip to next line dart 
end if 
end for 

count 5^ as a miss 
end while 



mentation in detail; in brief, we cast a line dart into 
the domain and intersect it with the disks of previ- 
ously accepted samples to generate a set of uncovered 
segments, then uniformly sample from those segments 
with a point dart. The user specifies an acceptable 
void volume F, the fraction of the domain uncovered 
by sample disks. We have a conservative stopping cri- 
teria based on the number of successive misses that 
usually achieves a smaller void volume. Figure 5 can 
be used as a guide for selecting V in four dimensions. 



4.1.1 Complexity 

The memory requirements are only 0{nd), which is 
what is required to represent the output point cloud 
because each sample has d coordinates. Only 2n -\- d 
floats are needed for scratch space for line segments g. 
We may generate and store only one fiat of one dart 
at a time. The runtime is 0{dn\ogn + nd^) per dart 
throw. The most significant feature is that the com- 
plexity does not suffer from a curse of dimensionality: 
there are no exponents containing d. The number of 
throws is a function of V and the miss rate; the com- 
munity appears to lack the tools to analytically bound 
the miss rate, but we show that for line darts it can 
be made to be reasonable, or at least more reasonable 
than the alternatives. 

Our approach is efficient because a line dart is more 
likely to intersect an uncovered region than a point 
dart. Only extremely simple and one-dimensional data 
structures are needed. The cost of throwing a line dart 
is nearly the same as a point dart. 



4.1.2 Output and Process 

The main drawback is that the output is not maximal, 
but its deviation can be estimated. A second potential 
drawback is that the process is not identical to MPS. 
There is no proof that the expected outputs are the 
same. Indeed, for a non-maximal sample, the probabil- 
ity of inserting the next sample point in a given disk- 
free subregion depends only on the subregion's area in 
MPS; but for our variant the probability depends also 
on how much prior disks cover the axis- aligned lines 
through the subregion. The main effect appears to be 
the order in which points are introduced, rather than 
their spatial distribution near maximality. Choosing 
the position of the next sample is dependent on many 
prior random decisions, so there are few noticeable pat- 
terns between one run and another with a different ran- 
dom number seed. 

MPS and Algorithm 1 follow a different process. The 
quality of the outcome, the point positions, does not 
appear to be sensitive to this, and we see little distinc- 
tion between the outputs using the standard measures 
of FFT spectrum, power, and anisotropy. Additionally, 
we are unaware of any formal definition of an ideal point 
cloud nor any proof that MPS produces it, so exactly 
matching the MPS process and its output are not strict 
requirements. The conventional goal in the graphics lit- 
erature is "blue noise," meaning no discernible axis or 
boundary aligned patterns, and something resembling 
a Heaviside step function for the mean radial power, as 
in Figure 10 right. 

4.1.3 Implementation Details 

k-d Tree To speed up the iteration over prior samples 
when generating segments g we use a, k-d tree. This 
allows us to prune samples whose disks are too far away 
to intersect the line of sj. A k-d tree requires little 
memory regardless of the number of dimensions. 

Line Segments g Beyond the k-d tree and storing 
the accepted samples, our only data structure is g, an 
array of segments in the line of sj. Each stored seg- 
ment is inside the domain but outside the sample disks 
processed so far. The set g may be implemented as a 
one-dimensional linked list storing the ith coordinate of 
the starts and ends of segments: g = [aoboaibi . . . aqbq] 
where q < n if the domain is convex. 

To update g for a new disk D{c) centered at sample 
c, we first compute g^ = D{c) H I = [6', a']. Using 
binary search in O(logn) time, we find the position in 
g where and a' should appear. If they are beyond ao 
or 6g, the disk does not intersect g. If b' and a' both 
lie between aj and 6j, the latter segment is split into 
two: ajb'a'bj. Otherwise, if only b' lies between aj and 
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6j, then a segment is trimmed by replacing bj by b'; 
a similar step applies for a^ Any endpoints between 
V and are covered by the disk and discarded. The 
updates are 0(1) time for a linked list. 

To choose a point from g uniformly, we sum the 
lengths of the segments L, choose a random number 
between and L, and use binary search to find the 
corresponding point on g. 



Maximality Estimates The user sets the remain- 
ing void volume V that is acceptable. Here we show 
how to translate that into a conservative stopping cri- 
teria. Denote the probability of hitting the void with 
a k-d dart by P^. Given a lower bound on P/^, we set 
m > [1/P/c]. We stop after m consecutive misses. 

We seek lower bounds on Pj^ in terms of V in order 
to find a sufficiently large value of m. For point darts, 
k = and Pq = V regardless of void shape. For higher 
dimensional darts, the worst shape for a void is a hy- 
percube with edge length b = Vd. By worst shape, we 
mean the shape that has the smallest Pk for a fixed 
assuming the entire domain is a hypercube. For a hy- 
percube void, the probability of hitting that void with 

d—k 

a k-d fiat is given by pk = . A k-d dart contains 
k = (f) k-d fiats and hence P^ = 1 - (1 -p/e)^'. This 
gives a lower bound on P/^, and a sufficient value of m. 
We may consider this as a MC estimation of the remain- 
ing void using a sample size of the last m + 1 k-d darts. 
As in any MC process, there is some variance, which 
decreases as k increases. Moreover, the remaining void 
is typically scattered throughout the domain; Tc is less 
than if the void volume was a single ball. 

4.1.4 Experimental Results 

Distribution Aspect Ratio The free radius, r/, is 
the disk radius, the minimum distance between any 
two samples. The coverage radius, Tc, is the maxi- 
mum distance between a domain point and its near- 
est sample. We deffne the distribution aspect ratio as 

= ^c/^/ ^ 1- This is a measure of maximality. 

To compute this for a point cloud, we used Qhull [1] 
to generate a Voronoi diagram. For each Voronoi vertex 
we retrieved the distance to its closest sample point: Tc 
is the maximum of these distances. 

Figure 5 shows the relation between the average dis- 
tribution aspect ratio and the acceptable void volume 
over four-dimensional samples of two different disk-free 
radii: rj =0.1 and rf = 0.05. Figure 7 shows run-times 
for generating the point clouds. Line darts consistently 
produced better results. 




Acceptable Remaining 



Figure 5: Achieved distribution aspect ratio for void 
volume threshold V in d = 4. Data points are averages 
over 25 experiments. Right shows more statistics for 
rf = 0.01, including standard deviation, sd. 



point darts and line darts. Figure 8 shows the number 
of points inserted over time. The expected void vol- 
ume V is related to the number of points; in practice 
the number of inserted points is a better indicator of 
maximality than our loose estimates of V based on suc- 
cessive misses. Line darts were able to generate larger 
samples for all V and d. 



Simple MPS ran out of memory (2GB) for r^< 0.025. 
For all cases < 1.3 r. . For Line darts, V = 0.01. 




0.04 0.05 0.06 0.07 0.08 
Disk-free Radius (r^) 



Figure 6: Time (blues) and memory (reds) for line darts 
compared to Simple MPS for the same acceptable void 
volume, V=le-2. 



Efficiency by Method Figure 6 compares the per- 
formance of traditional point darts and line darts (this 
work) to Simple MPS [5]. Recah Simple MPS is based 
on a ffat quadtree, and is currently the fastest and most 
memory-efficient of the provably- correct MPS methods. 
Our method is attractive at large values of rj for its 
speed, and at low values of r/ for its memory consump- 
tion. For example, for r/ < 0.025, we generated 4M 
points in half an hour using 107 MB of memory, while 
Simple MPS ran out of memory at 2 GB. 



Speed of Approaching Maximality We tested our 
code over 2, 4, 10, and 30-dimensional domains using 



Output Quality We measure the quality of the dis- 
tribution of 2-d output points using the PSA spectrum 
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E lE+2 




Distribution Aspect Ratio 



Figure 7: V threshold effects on time (blues) and sam- 
ple size (greens) for line and point darts. Data points 
are labeled with log^^g ^ decreases, the sample 

approaches maximality and distribution aspect ratio 1. 
Simple MPS [5] dashed lines are for a maximal distri- 
bution. 



analysis tool [27]. Using point darts, our process is the 
same as classic dart-throwing, so we it for our stan- 
dard of correct output. Figures 9 and 10 compare the 
outputs' blue-noise properties; the difference between 
point and line darts was insignificant, at least for these 
three metrics. 



4.2 Depth of Field With Antiahasing 

k-d darts can be used for fast and high-quality ren- 
dering of depth-of-field (DOF) effects in computer- 
synthesized images. Mathematically, computing a 
pixel's color in the presence of DOF can be expressed as 
a four-dimensional integral over the pixel's spatial (x, y) 
and lens aperture (u^v) dimensions. In most high- 
quality renderers, this is calculated using Monte Carlo 
integration over many point samples. This method suf- 
fers from a low rate of convergence; reducing noise for a 
good quality image usually requires a very large num- 
ber of samples per pixel, k-d darts offer the promise 
of faster convergence with low noise. Instead of us- 
ing point samples for reconstruction, we use 1-d (line) 
darts, thrown in the 4-d (x, v) space. 

We use Latin Hypercube Sampling (LHS) or jittered 
sampling for each dimension [2] . Given that our sample 
space is four-dimensional, each line dart consists of four 
line flats. We select n points, i.e. An flats. Each line 
requires a fixed location in 3-d and a variable fourth 
dimension. 

We compute coverage for these darts using a method 
inspired by Gribel et al.'s work [10, 9] on rendering 
motion blur. Gribel et al. fix x and y and shoot line 
samples in the time domain; instead, we use line darts 
that shoot line flats in all spatial dimensions to com- 
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Figure 8: Line darts approach maximality faster than 
point darts, as measured by the number of inserted 
points in a given run time. 




(a) Line darts 




(b) Point darts 



Figure 9: FFT spectra (right) for the relaxed MPS 
point clouds (left) generated by line darts (top) and 
point darts (bottom). 



pute depth-of-field blur, 
dimensional problem. 



We also sample a higher- 
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(a) Line darts 



• Tzeng et al. use non-random dart locations, but 
we randomly position darts. 

• Tzeng et al.'s work specifically targets the GPU 
pipeline; we do not discuss (or consider) imple- 
mentation details. 




frequency 



180 360 540 720 

frequency 



(b) Point darts 



Figure 10: Radial anisotropy and mean power esti- 
mates for line darts (top) and point darts (bottom), 
averaged over ten samplings. 

4.2.1 Contrast to Tzeng et al. 

Tzeng et al. [30] considered line sampling for DOF. 
They only consider the two {u, v) dimensions, not the 
four (u^v^x^y) dimensions as we do. Further, they take 
advantage of the structure that the (ii, v) subspace of 
interest is uniformly circular. This reduces 2-d (u^v) 
sampling to 1-d sampling of lines through the origin, 
pinwheel sampling by angle. In their implementation, 
for each pixel, they fix x and y and use a pinwheel 
of line samples to vary u and v. The resulting DOF 
had high performance when compared to point sam- 
pling strategies, with low noise. However, due to both 
the pinwheel configuration and the fixing of x and y for 
each pixel's sample, strobing artifacts tended to occur 
in regions with high frequency changes. The formula- 
tion we present in this paper differs in the following 
ways: 

• Both implementations address DOF, but we also 
address antialiasing. 

• Tzeng et al.'s line darts are all radial and speci- 
fied completely by their angle; they live in 1-d {0) 
space. In this paper, we use axis-aligned darts in 
4-d (x, y, v) space. The positive and negative 
consequences follow: 

— Tzeng et al. exhibits screen-space aliasing be- 
cause X and y are fixed. 

— Tzeng et al. achieves higher quality DOF for 
the same number of samples, because 1-d 
spaces require fewer samples to cover than 4- 
d spaces. 



4.2.2 Triangle Edge Equations in 4-d 

We now describe the triangle equations and how we 
create a line sample. For a given triangle, we start 
by computing a signed radius of Circle of Confusion 
(CoC) for each vertex, obtained using the following ex- 
pression [13]: 



CoC = A 



where A and / are the camera aperture and focal 
length, respectively, and Zf and z indicate the respec- 
tive depths of the focal plane and the given vertex. 
Note z is simply the w coordinate of the vertex in clip 
space. Now that we have the circle of confusion for each 
triangle's vertex, we can begin formulating a sampling 
strategy for each pixel. 

Given a set of coordinates on the lens u and we as- 
sume a linear apparent motion of the vertex screen co- 
ordinates. For a given screen space vertex i G {0, 1, 2}, 
with coordinates {xi^yi)^ circle of confusion q, and 
u^v ^ [—0.5,0.5], we have 



Xj^ — Xi ~\~ CiU 

Vi =yi + CiV 



(2) 
(3) 



For each pixel we have a four- dimensional space 
(x^y^u^v), where x and y are the subpixel regions, 
and u and v are coordinates on the lens. To get a 
realistic and noise-free image, we seek to sample this 
space uniformly in an effective manner. Since we have 
four varying dimensions, we choose to use four different 
hypercubes for our Latin Hypercube Sampling (LHS). 
Each hypercube chooses three dimensions within which 
to sample, and one that will vary in our next stage. 

For each triangle, we can consider the edge equations 
in this four-dimensional space to be the following: we 
substitute Equations (2) and (3) into the equations for 
testing whether a point (x, y) lies within one of the 
triangle edges i. Let ESi represent the edge-sum at 
point {x^y) for the ith edge, between vertices i and j, 

ESiix, y, u, v) = iy- yl){x-, - x^) - (x - x^M ' vD- 
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I Given yi and Vi 
I >^ 
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E(u),= P - Q(,u) 







(a) Two line darts in (x, u) (b) Resolving occlusion depth 

Figure 11: Our technique for computing analytical cov- 
erage using line darts. In this example, consider two 
possible line flats in the x and u directions. The do- 
main can be transformed and we can test the triangles 
for occlusion along the line flat as shown in (a). Depth 
resolution per flat is shown in (b). 



Expanding, we get equivalent right-hand sides 

^ {y-Vi- Q^) {xj -Xi^ u{cj - Ci)) 

- {x-Xi- Ciu) {yj -Vi^ v{cj - Ci)) 
^ ESi{x,y, 0,0) 

+ u{cj{y-yi) -Ci{y-yj)) 

- V {Cj{x - Xi) - Ci{x - Xj)) . 

This simplifies to 

ESi{x, y, u, v) = Ci{x, y) + uAi{y) - vBi{x) > 0. 

Here we have four dimensions of variability 
(x^y^u^v). In conventional point sampling, one would 
randomly (or pseudo-randomly) select a set of points in 
this space. We instead decide to employ line sampling 
using our Latin hypercubes. To illustrate this concept, 
consider fixing three of these dimensions with points on 
our hypercube from LHS: select xi^yi and Ui for x,y 
and u. Now we have a line equation as follows: 

ESi{v) = Ci{xi,yi) ^uiAi{yi) - vBi{xi), 

which simplifies to 

ESi{v) = P-vQ. 

This line equation is easy to analytically solve for 
each of our edge equations and determines where line 
coverage exists. The same concept can be extended to 
our other hypercube setups, fixing three of the dimen- 
sions and varying the last. 

4.2.3 Analytical Coverage in the Hypercube 
Domain 

Knowing the edge equations in our 4-d domain, we can 
now compute their coverage along a line sample. Fig- 



Table 1: Performance of our point vs. line darts. 
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1024 points 
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256 points 


1024 points 


16 lines 


30 lines 


time (s) 


214 


853 


157 


292 



ure 11 summarizes our technique for determining cov- 
erage. 

We instantiate a set of line flats along each of our four 
hypercube configurations. A line dart is the combina- 
tion of four different line flats (one in the direction of 
each of the four dimensions of the domain) with initial 
points chosen using our LHS. 

Rendering consists of testing each incoming triangle 
against potentially covered line flats. For each pixel 
sample in the triangle's bounding box, we use equations 
from Section 4.2.2 to transform triangle edges to test 
for the correct hypercube domain. 

To finish the calculation we follow Gribel et al.'s ap- 
proach [10] to construct and resolve line darts. 

For each line flat, we analytically compute its seg- 
ment covered by the triangle. A per-line-sample queue 
stores the color and depth of covered segments. Once 
all triangles have been processed, we resolve the final 
color for each sample. We sweep across the flat while 
aggregating triangles closest in depth, and then use all 
pixel samples to compute the final color for the pixel. 

4.2.4 Implementation and Results 

To test our formulation, we built a simple CPU-based 
renderer. It is capable of rendering scenes using tradi- 
tional triangle rasterization. We integrated two addi- 
tional capabilities: 

• A stochastic sampler based on point darts. 

• A sampler based on line darts. 

The two noise artifacts that typically occur using a 
DOF scheme are noise in blurry regions and noisy alias- 
ing artifacts near the point in focus. For scenes far 
away from the focal plane, line flats in the x and y di- 
rection become particularly narrow when compared to 
the sampling space, and this causes additional noise. 
In these regions the length of flats in the x and y di- 
mensions are significantly smaller than flats in either 
the u or V direction. 

If we consider all such flats to have equal contribu- 
tions, this results in the x and y samples adding a no- 
ticeable amount of noise into our system in blurry re- 
gions. To address this, rather than considering all line 
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fiats to have equal contribution, we select a weight con- 
stant (a = 0.2 in our case) such that x and y line flats 
are scaled by a and contribute less to the scene. 

This weight can be modified based on the aperture 
of the lens for the scene. In cases where the aperture 
is small, contributions from x and y are deemed more 
important (for antialiasing), and the weight is adjusted 
accordingly. A more accurate dart sampling method 
might also consider the ratios of flat lengths per tile 
(worst case) and decide on an a weight accordingly. 
Research is needed to determine the optimal amount 
for each fiat to contribute. Still, our simple heuristic 
seems to be effective. 

Figure 12 compares two scenes rendered with our 
point dart and line dart techniques. Renderings based 
on k-d darts are virtually free of noise and aliasing ar- 
tifacts. Point darts, however, retain noticeable aliasing 
artifacts even with 1024 darts. Although 16 line darts 
produce a bit more noise in unfocused regions than 256 
point darts, 30 line darts has better quality than 256 
point darts and around the same quality as 1024 point 
darts in unfocused regions, and no noticeable aliasing 
artifacts in focused areas. 

Table 1 shows the performance of our samplers. 
Clearly, throwing one line dart is more expensive than 
throwing one point dart. However, fewer are needed, 
and correctly weighted line darts converge more quickly 
to a less noisy image without aliasing artifacts. 




2 "1 

(a) Circular Parabola 




2 -1 

(b) Planar Cross 



Figure 13: Plots of the two surrogate model test func- 
tions. 
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4.3 Probability of Failure 

Uncertainty quantification usually explores a vast high- 
dimensional space with a limited budget of sample 
points. Efficiency is crucial because typically the func- 
tion evaluation is expensive and more sample points 
are desired than we can afford. Surrogate models cre- 
ate a cheaper response surface that is evaluated instead. 
Even for very cheap response surfaces, if the failure re- 
gion is small enough Monte Carlo (MC) sampling will 
not estimate it accurately. Here we show that k-d darts 
can improve MC efficiency. 

We test the "circular parabola" (Equation (4)) and 
"planar cross" (Equation (5)) surrogate models; see 
Figure 13. 

d 

y{x) = ^{2x,-lf, 0<x,<l. (4) 
i=i 

l/d 



y{x) = 



n 



(1 + cos(27rXi)) 



, <Xi <1. (5) 



Failure is defined as the function value below some 
constant threshold: y{x) < yt. The shape of the fail- 
ure region is different for the two test functions: a disk 



Figure 14: Speedup of line darts over point sampling 
for estimating the probability of failure of two analytic 
response surfaces. 



for the parabola and a fattened plus-sign for the planar 
cross; see Figure 13 and also Figure 1. For uniform dis- 
tributions, the probability of failure is the fraction of 
the domain volume where y{x) < yt. We choose yt so 
the probability of failure is exactly 10~^ or 10~^. We 
estimate the failure volume using line darts. For a line 
flat, we find the roots of the single variable equation 
y{xi) = yt- The length of the line segment between the 
two roots (if real roots exist) is used to estimate the 
volume of the failure region. Figure 14 demonstrates 
the benefit of line darts over conventional point sam- 
pling in reducing the time required to achieve a given 
accuracy level. The complexity of root-finding reduces 
the performance of line darts. However, for both test 
functions, line darts were better than point darts in 
dimensions up to 15. 
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(a) Point darts, 256 per pixel. 




(b) Point darts, 1024 per pixel. 




(c) Line darts, 16 per pixel. 




(d) Line darts, 30 per pixel. 



Figure 12: Depth Of Field (DOF) images using conventional point sampling (a-b) vs. our k-d darts (c-d). k-d darts 
produce high-quality, antialiased images. The third column from the left shows a close up of an extremely blurry 
region, the Cessna's tail. In regions close to the focal plane we see some aliasing artifacts for point darts but not for 
line darts, e.g. the Cessna's prop cone, the junction of its body and wing, and the transition shades on the green 
teapot spout. Furthermore, line darts tend to be faster for the same quality blur; see Table 1. 



5 Accuracy Experiments 

5.1 Problem Motivation 

We provide some experimental results on the accu- 
racy of darts for the canonical Monte Carlo problem 
of estimating the volume of an object in high dimen- 
sions. (Volume estimation is in the same category as 
the probability-of-failure problem in Section 4.3.) In 
particular, we seek to show that the method produces 
good estimates, regardless of the size, shape, dimen- 
sion and orientation of the object, and regardless of 



the dimension and orientation of the darts. The aver- 
age estimate should be close to the true estimate, and 
the higher moments of the estimates should be low. We 
design our experiments to show the effects (if any) of 
the following factors. 

• the dimension of the object. 

• /c, the dimension of the dart. Of particular inter- 
est is comparing our results to standard MC point 
sampling, /c = 0. 

• 5, the squish factor of the object, which controls 
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its aspect ratio. 

• r, the number of rotations of the object. This 
aUows us to compare axis-ahgned objects to un- 
ahgned ones. 

• Axis-ahgned darts vs. unahgned darts. 

We perform experiments of n flats over the prior pa- 
rameters, as described in Table 2. Note that we keep 
the number of flats constant, rather than the number 
of darts, because the computational expense is more 
closely tied to the number of flats for objects with an 
analytic expression, and because middle-dimensional 
darts have many more flats than high and low dimen- 
sional ones. 



Table 2: k-d dart parameter study 
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The darts are aligned with the coordinate axis except 
for some d = 2 experiments. We repeat each parameter 
combination 1000 times, N = 1000. Squish parameters 
s are symmetric around 1 with respect to ellipse 
volume: e.g. for (i = 2, 5 = 1/2 and s = 2 define 
ellipses with the same volume. 



5.2 Object Generation 

Instead of a spherical object, we estimate the volume 
of an ellipse (a.k.a. ellipsoid), randomly oriented and 
squished. An ellipse provides enough generality to test 
the factors in Section 5.1, but enough simplicity to iso- 
late numerical from methodology issues. In particular, 
we choose an ellipse because it is possible to analyt- 
ically calculate the volume of an ellipse's intersection 
with a k-d dart. 

Our object is a (i- dimensional ellipse centered at the 
origin. We construct it as follows. We start with a 
(i-sphere centered at the origin with radius 1. This 
fits in an origin-centered cube with side length 2, the 
two-cube. We ensure the final ellipse also lies in the 
two-cube. 

• s: squish factor. We scale the ellipse along the x- 
axis by multiplying its x-extent by s. The sphere 
has 5 = 1. Note 5 < 1 gives thin, coin-shaped 
objects. For 5 > 1, we then shrink the ellipse so 
it fits in the two-cube: multiply all coordinates by 



a factor of 1/s. The net effect is keeping the x- 
coordinate fixed and scaling the other axes by 1/s. 
This gives needle-shaped objects. 

• r: number of rotations. We perform r random ro- 
tations in sequence. Each is a Givens rotation with 
a random pair of coordinate indices i and j, and a 
random angle G [0,7r]: multiply coordinates by 
the identity matrix with the 2x2, {i, j} submatrix 
[1 0; 1] replaced by [cos^ -sin^; sin^ cos^]. 

5.3 Dart Generation 

Most implementers will choose axis-aligned darts for 
three reasons. First, it is easy to distribute aligned 
darts uniformly, which ensures that the expected mean 
of the function estimates is accurate. Second, it is easi- 
est to implement aligned darts, since it involves simply 
fixing coordinate values. Third, in many cases it is most 
efficient because we may obtain an expression for the 
underlying function along a dart by substituting in the 
fixed coordinate values. However, for completeness we 
provide some experimental results on the accuracy of 
unaligned darts. We shoot k-d darts into the two-cube 
as follows: 

• A point dart {k = 0) is generated by selecting a 
random point. Each of the d coordinates is chosen 
independently, and uniformly in [0, 1]. 

• Aligned darts are generated by their fiats. Each 
flat has a unique combination of d — k fixed co- 
ordinate indices; the remaining k coordinates are 
allowed to vary. The coordinates for the fixed in- 
dices are chosen independently and uniformly as 
for point darts. 

• Unaligned flats are generated so that the orienta- 
tion of the fiats is uniformly random. The only 
experimental setting where we generate unaligned 
flats is for A: = 1 and d = 2, line darts in the plane. 
We choose angle G [0,7r], which determines the 
orientation of the flat. Any line that intersects 
the square crosses one of its main diagonals. (It 
is guaranteed to cross the diagonal to which it is 
more perpendicular, which depends only on ^.) We 
pick a point p uniformly at random along the ap- 
propriate diagonal. We now have a point and an 
angle, which defines a line flat. For random darts, 
the second flat is a line perpendicular to the first 
line, passing through some random point q of the 
other diagonal. 

Aligned 1-d darts are labeled "k=la," random flats 
are labeled "k=lr," and random darts, pairs of orthog- 
onal flats, are labeled "k=lo" in the top two rows of 
Figures 15 and 16. 
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5.4 Object-Dart Intersection 

For point darts, the volume estimation is the fraction 
of darts that landed inside the ellipse, multiplied by 
the volume of the two-cube sampling domain, 2^. For 
k > darts, instead of this discrete ratio we average 
the geometric fraction of each dart inside the ellipse 
object. The details of these calculations follow. 

For point darts, we simply back-project the points 
to the domain of the sphere: apply the inverse Givens 
rotations to the dart's point in reverse order; then scale 
the x-coordinate by I/5 (or, for 5 > 1, all the other 
coordinates by s). If the distance from the transformed 
dart to the origin is less than 1 it is inside the sphere, 
and the original dart is inside the ellipse. 

For k-d darts, we back-project their hyperplanes into 
the sphere domain, where we can calculate the volume 
of intersection analytically, and then forward-weight it 
by the scaling. 

To back-project a fiat, we back-project A: + 1 points 
spanning the fiat. Each dart has d — k fixed coordinates 
and k free coordinates. We pick spanning point po with 
for all of its free coordinates, and spanning point p^>o 
with 1 for its ith free coordinate and for its other 
free coordinates. Each pi is back-projected to qi using 
the same procedure as for a point dart. The k vectors 
from qo to {qi>o} span the transformed fiat, but are no 
longer orthonormal because of the final scaling step, 
so we must reconstruct an orthonormal basis. Now we 
are ready to calculate volumes, using forward trans- 
formations. We calculate the distance from the flat to 
the origin. This tells us the radius of the /c- dimensional 
subsphere that is the intersection of the fiat with the d- 
sphere. We compute the volume of this subsphere. We 
multiply this volume by the sum of the x-components 
of the orthonormal basis, which gives the volume of 
the (unrotated) ellipse of intersection. The (forward) 
rotations do not affect the volume so are skipped. 

For unaligned line-darts in the plane, the distance 
from the origin is easy to measure. We use a process 
similar to the prior paragraph but it is a little easier 
because the 1-sphere is simply a line segment. 

5.5 Results 

We plot the mean of the absolute value of the relative 
error, |mean — true | /true, vs. the number of flats n 
in Figure 15. We plot the histograms of the ratios of 
estimated / true volume for n = 10^ in Figure 16. Each 
subfigure shows results for all k for some combination 
of the other parameters. We observe the following. 

In Figure 15 the experimental slopes for all k and 
d are about —1/2, in agreement with theory (the n 
to the power —1/2 in Equation (1)). The accuracy is 
insensitive to the orientation of the object. 



In Figure 16 the histograms are all sharply peaked 
at the true value. This shows that the variations, the 
higher order moments of the estimates, are reasonable. 

The major trend of these figures is that the accuracy 
of the estimates improves with k. Moreover, the larger 
the A:, the smaller the variation of the estimate and the 
sharper the peak near the true value. 

For small- volume objects, aligned darts are more ac- 
curate than unaligned darts. This is illustrated by the 
red curves in the top-right and top-left subfigures in 
Figures 15 and 16. We were initially surprised by this, 
but the explanation is that unaligned fiats are shorter 
on average than aligned ones, because they might clip a 
corner of the square rather than being the side length of 
the square, so they more often miss the object. (They 
will be even shorter as the dimension of the space in- 
creases.) For moderate-volume objects, the accuracy 
is about the same regardless of object orientation. For 
unaligned fiats, it appears that our n was large enough 
that using pairs of orthogonal fiats or independently- 
random individual fiats does not make much difference. 
Our conclusion is that aligned darts are universally bet- 
ter when the sample domain is a square. 

The accuracy is primarily sensitive to the volume of 
the object and, secondarily, to the squish value. Higher 
dimensional darts are better than lower dimensional 
ones, and the advantage is more pronounced for small- 
volume objects. This can been seen by considering the 
second- from-bottom row in Figures 15. To see the vol- 
ume dependence, note that the lines are closer together 
for 5 = 1, and farther apart for larger and smaller 
squish factors. In low dimensions, 2 and 3, the smaller 
the volume of the object, the less accurate are all es- 
timates, for all dimensions k. Mo derate- dimensional 
darts in high dimensions, e.g. d = 10 and /c = 4, also 
exhibit this trend. However, 9-d darts in 10-d space 
have the same accuracy regardless of the volume or 
squish, because they are close to the dimension of the 
underlying space and they more fully span it. For ex- 
ample, a dart with k = d always evaluates to the true 
value. That is, the advantage of higher- dimensional 
darts over lower-dimensional darts is greater for small 
objects in high dimensions, which is where we are ad- 
vocating their use. 

A secondary phenomena is that the estimate is 
slightly more accurate for squish s = 1/10 than for 
s = despite the volume being the same. This 

can be seen by careful examination of the height of the 
lines in the d = 10, varying squish row in Figure 15. 
For instance, the mean error is 10% lower in the first 
column than the fifth for n = 10^ and k = 4. This sup- 
ports our intuition that darts are effective at hitting 
thin, coin-shaped regions. Since ^ 1-3, a 10-d 

object with this squish factor is actually roundish and 
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d = 2, varying squish 





Figure 15: Mean error of volume estimation, |mean — true | /true by n. See also Table 2 and Figure 16. 



not very needle-shaped. Preliminary experiments show 
that the accuracy gained by increasing /c is a compli- 
cated function of s and the volume. High dimensional 
darts have more advantage over low dimensional ones 
for very small and sharp needle-shaped objects, com- 
pared to their advantage for coin-shaped objects. 

6 Conclusions 

In this work we introduce k-d darts as a conception of 
higher-dimensional sampling. We described a k-d dart 



framework for general dimension /c, and then demon- 
strated efficiency and accuracy over three applications 
using using /c = 1, and accuracy for one application us- 
ing k > 1. In particular, we showed that darts produce 
accurate estimates of the volume of an object regard- 
less of the dimension, orientation and aspect ratio of 
the object. Axis- aligned darts are universally prefer- 
able to unaligned ones for sampling square domains, 
and we expect this to extend to hyper-rectangles, e.g. 
bounding boxes. 

In principle, high- aspect-ratio (small- volume) ob- 
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Figure 16: Volume estimation histograms, estimate / true by frequency, for n = 10^. "la" is axis-aligned darts, "Ir" 
is random-angle lines, and "lo" is random-orientation darts, pairs of orthogonal flats. See also Table 2 and Figure 15. 



jects are more efficiently sampled by higher- /c darts. 
Demonstrating that efficiency for applications requires 
either analytical expressions, as in the toy sphere- 
volume problem in the Introduction; or efficient numer- 
ical techniques for evaluating the underlying function 
along higher dimensional flats. In future work we plan 
to explore a numerical technique using recursive sam- 
pling designs by dimension. 

For maximal Poisson-disk sampling, line darts are 
helpful in getting close to maximality in high dimen- 



sions. Below a given acceptable void ratio, they are 
more efficient than point darts. In terms of bounding 
the distance from a domain point to the nearest sam- 
ple point, we are actually closer to maximality as the 
dimension increases. Any difference in the distribution 
of points produced by classical maximal Poisson-disk 
sampling and our line darts is small by the standard 
measures. Our line dart algorithm is efficient with re- 
spect to memory usage, which enables the production 
of larger samples in higher dimensions. 
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For depth-of- field, using generalized k-d darts gives 
us a high-quality noise-free image without aliasing. Al- 
though each 1-d dart requires more processing than 
a point sample, we only need a few of them to ren- 
der a high-quality image. This results in our k-d dart 
method outperforming point sampling. We suggest 
sampling over a higher-dimensional space to render 
both depth-of-field and motion blur in animations. We 
suggest exploring weighted 1-d darts, where scene in- 
formation determines which fiats contribute more to a 
scene. 

For uncertainty quantification, k-d dart Monte Carlo 
sampling can be more efficient than point sampling. 
The key for all these applications is exploiting the prob- 
lem structure to take advantage of what k-d darts pro- 
vide. 
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